grosbeak investigate Z chromosome

0.1 Investigate genotype patterns in hybrids

Code
library(triangulaR)
Loading required package: ggplot2
This is triangulaR v.0.0.1

          /\
         /  \
        /    \
       /______\

Usage information is available at: https://github.com/omys-omics/triangulaR/ 

Please cite the following if you use triangulaR in a publication: 
Code
library(vcfR)

   *****       ***   vcfR   ***       *****
   This is vcfR 1.14.0 
     browseVignettes('vcfR') # Documentation
     citation('vcfR') # Citation
   *****       *****      *****       *****
Code
library(ggplot2)
library(reshape2)
library(adegenet)
Loading required package: ade4

Attaching package: 'ade4'
The following object is masked from 'package:triangulaR':

    triangle.plot

   /// adegenet 2.1.10 is loaded ////////////

   > overview: '?adegenet'
   > tutorials/doc/questions: 'adegenetWeb()' 
   > bug reports/feature requests: adegenetIssues()
Code
library(StAMPP)
Loading required package: pegas
Loading required package: ape
Registered S3 method overwritten by 'pegas':
  method      from
  print.amova ade4

Attaching package: 'pegas'
The following object is masked from 'package:ape':

    mst
The following object is masked from 'package:ade4':

    amova
The following objects are masked from 'package:vcfR':

    getINFO, write.vcf
Code
#read in data
v<-read.vcfR("~/Desktop/grosbeak.rad/grosbeak.filtered.snps.vcf.gz")
Scanning file to determine attributes.
File attributes:
  meta lines: 44796
  header_line: 44797
  variant count: 51595
  column count: 147

Meta line 1000 read in.
Meta line 2000 read in.
Meta line 3000 read in.
Meta line 4000 read in.
Meta line 5000 read in.
Meta line 6000 read in.
Meta line 7000 read in.
Meta line 8000 read in.
Meta line 9000 read in.
Meta line 10000 read in.
Meta line 11000 read in.
Meta line 12000 read in.
Meta line 13000 read in.
Meta line 14000 read in.
Meta line 15000 read in.
Meta line 16000 read in.
Meta line 17000 read in.
Meta line 18000 read in.
Meta line 19000 read in.
Meta line 20000 read in.
Meta line 21000 read in.
Meta line 22000 read in.
Meta line 23000 read in.
Meta line 24000 read in.
Meta line 25000 read in.
Meta line 26000 read in.
Meta line 27000 read in.
Meta line 28000 read in.
Meta line 29000 read in.
Meta line 30000 read in.
Meta line 31000 read in.
Meta line 32000 read in.
Meta line 33000 read in.
Meta line 34000 read in.
Meta line 35000 read in.
Meta line 36000 read in.
Meta line 37000 read in.
Meta line 38000 read in.
Meta line 39000 read in.
Meta line 40000 read in.
Meta line 41000 read in.
Meta line 42000 read in.
Meta line 43000 read in.
Meta line 44000 read in.
Meta line 44796 read in.
All meta lines processed.
gt matrix initialized.
Character matrix gt created.
  Character matrix gt rows: 51595
  Character matrix gt cols: 147
  skip: 0
  nrows: 51595
  row_num: 0

Processed variant 1000
Processed variant 2000
Processed variant 3000
Processed variant 4000
Processed variant 5000
Processed variant 6000
Processed variant 7000
Processed variant 8000
Processed variant 9000
Processed variant 10000
Processed variant 11000
Processed variant 12000
Processed variant 13000
Processed variant 14000
Processed variant 15000
Processed variant 16000
Processed variant 17000
Processed variant 18000
Processed variant 19000
Processed variant 20000
Processed variant 21000
Processed variant 22000
Processed variant 23000
Processed variant 24000
Processed variant 25000
Processed variant 26000
Processed variant 27000
Processed variant 28000
Processed variant 29000
Processed variant 30000
Processed variant 31000
Processed variant 32000
Processed variant 33000
Processed variant 34000
Processed variant 35000
Processed variant 36000
Processed variant 37000
Processed variant 38000
Processed variant 39000
Processed variant 40000
Processed variant 41000
Processed variant 42000
Processed variant 43000
Processed variant 44000
Processed variant 45000
Processed variant 46000
Processed variant 47000
Processed variant 48000
Processed variant 49000
Processed variant 50000
Processed variant 51000
Processed variant: 51595
All variants processed
Code
v
***** Object of Class vcfR *****
138 samples
1027 CHROMs
51,595 variants
Object size: 532.6 Mb
2.987 percent missing data
*****        *****         *****
Code
colnames(v@gt)
  [1] "FORMAT"                 "P_hybrid_44696"         "P_hybrid_44703"        
  [4] "P_hybrid_44707"         "P_hybrid_44708"         "P_hybrid_44709"        
  [7] "P_hybrid_44712"         "P_hybrid_44762"         "P_hybrid_44771"        
 [10] "P_hybrid_44781"         "P_hybrid_45171"         "P_hybrid_45173"        
 [13] "P_hybrid_45174"         "P_ludovicianus_11998"   "P_ludovicianus_21721"  
 [16] "P_ludovicianus_25286"   "P_ludovicianus_26595"   "P_ludovicianus_33988"  
 [19] "P_ludovicianus_34776"   "P_ludovicianus_34779"   "P_ludovicianus_34782"  
 [22] "P_ludovicianus_34830"   "P_ludovicianus_44704"   "P_ludovicianus_44705"  
 [25] "P_ludovicianus_44706"   "P_ludovicianus_44710"   "P_ludovicianus_44711"  
 [28] "P_ludovicianus_44713"   "P_ludovicianus_44714"   "P_ludovicianus_44715"  
 [31] "P_ludovicianus_44716"   "P_ludovicianus_44719"   "P_ludovicianus_44720"  
 [34] "P_ludovicianus_44722"   "P_ludovicianus_44723"   "P_ludovicianus_44726"  
 [37] "P_ludovicianus_44727"   "P_ludovicianus_44729"   "P_ludovicianus_44730"  
 [40] "P_ludovicianus_44731"   "P_ludovicianus_44732"   "P_ludovicianus_44733"  
 [43] "P_ludovicianus_44734"   "P_ludovicianus_44737"   "P_ludovicianus_44738"  
 [46] "P_ludovicianus_44739"   "P_ludovicianus_44740"   "P_ludovicianus_44741"  
 [49] "P_ludovicianus_44742"   "P_ludovicianus_44743"   "P_ludovicianus_44744"  
 [52] "P_ludovicianus_44745"   "P_ludovicianus_44746"   "P_ludovicianus_44747"  
 [55] "P_ludovicianus_44748"   "P_ludovicianus_44749"   "P_ludovicianus_44753"  
 [58] "P_ludovicianus_44754"   "P_ludovicianus_44761"   "P_ludovicianus_44775"  
 [61] "P_melanocephalus_34890" "P_melanocephalus_43110" "P_melanocephalus_43276"
 [64] "P_melanocephalus_44346" "P_melanocephalus_44347" "P_melanocephalus_44471"
 [67] "P_melanocephalus_44651" "P_melanocephalus_44660" "P_melanocephalus_44666"
 [70] "P_melanocephalus_44676" "P_melanocephalus_44677" "P_melanocephalus_44678"
 [73] "P_melanocephalus_44679" "P_melanocephalus_44680" "P_melanocephalus_44681"
 [76] "P_melanocephalus_44683" "P_melanocephalus_44684" "P_melanocephalus_44685"
 [79] "P_melanocephalus_44686" "P_melanocephalus_44687" "P_melanocephalus_44688"
 [82] "P_melanocephalus_44689" "P_melanocephalus_44692" "P_melanocephalus_44693"
 [85] "P_melanocephalus_44694" "P_melanocephalus_44695" "P_melanocephalus_44697"
 [88] "P_melanocephalus_44699" "P_melanocephalus_44752" "P_melanocephalus_44760"
 [91] "P_melanocephalus_44763" "P_melanocephalus_44764" "P_melanocephalus_44765"
 [94] "P_melanocephalus_44766" "P_melanocephalus_44769" "P_melanocephalus_44770"
 [97] "P_melanocephalus_44772" "P_melanocephalus_44773" "P_melanocephalus_44774"
[100] "P_melanocephalus_44776" "P_melanocephalus_44777" "P_melanocephalus_44778"
[103] "P_melanocephalus_44779" "P_melanocephalus_44780" "P_melanocephalus_44782"
[106] "P_melanocephalus_44787" "P_melanocephalus_44788" "P_melanocephalus_44789"
[109] "P_melanocephalus_44790" "P_melanocephalus_44791" "P_melanocephalus_44792"
[112] "P_melanocephalus_44793" "P_melanocephalus_44794" "P_melanocephalus_44795"
[115] "P_melanocephalus_44798" "P_melanocephalus_44799" "P_melanocephalus_44801"
[118] "P_melanocephalus_44802" "P_melanocephalus_44803" "P_melanocephalus_44804"
[121] "P_melanocephalus_44805" "P_melanocephalus_44806" "P_melanocephalus_44808"
[124] "P_melanocephalus_44809" "P_melanocephalus_44810" "P_melanocephalus_44837"
[127] "P_melanocephalus_44840" "P_melanocephalus_44843" "P_melanocephalus_44845"
[130] "P_melanocephalus_44846" "P_melanocephalus_44847" "P_melanocephalus_44848"
[133] "P_melanocephalus_44853" "P_melanocephalus_44854" "P_melanocephalus_45175"
[136] "P_melanocephalus_45200" "P_melanocephalus_45232" "P_melanocephalus_45709"
[139] "P_melanocephalus_45926"
Code
#read in sampling dataframe
samps<-read.csv("~/Desktop/grosbeak.data.csv")
samps<-samps[samps$passed.genomic.filtering == "TRUE",]
#make popmap
pm<-data.frame(id=samps$sample_id,pop=c(rep("P1", times=9),rep("P0",times=8),rep("hyb",times=121)))

# Create a new vcfR object composed only of sites above the given allele frequency difference threshold
fixed.vcf <- alleleFreqDiff(vcfR = v, pm = pm, p1 = "P0", p2 = "P1", difference = 1)
[1] "266 sites passed allele frequency difference threshold"
Code
# Calculate hybrid index and heterozygosity for each sample. Values are returned in a data.frame
hi.het.test <- hybridIndex(vcfR = fixed.vcf, pm = pm, p1 = "P0", p2 = "P1")
[1] "calculating hybrid indices and heterozygosities based on 266 sites"
Code
# View triangle plot
triangulaR::triangle.plot(hi.het.test)
Warning in geom_segment(aes(x = 0.5, xend = 1, y = 1, yend = 0), color = "black"): All aesthetics have length 1, but the data has 138 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Warning in geom_segment(aes(x = 0, xend = 0.5, y = 0, yend = 1), color = "black"): All aesthetics have length 1, but the data has 138 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Warning in geom_segment(aes(x = 0, xend = 1, y = 0, yend = 0), color = "black"): All aesthetics have length 1, but the data has 138 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.

Code
#see where the fixed differences are coming from:
table(fixed.vcf@fix[,1])

VZSJ01000159.1 VZSJ01000241.1 VZSJ01000254.1 VZSJ01000265.1 VZSJ01000270.1 
             3              2              3              5             11 
VZSJ01000272.1 VZSJ01000279.1 VZSJ01000280.1 VZSJ01000283.1 VZSJ01000285.1 
             3              9              3              4              3 
VZSJ01000299.1 VZSJ01000301.1 VZSJ01000305.1 VZSJ01000337.1 VZSJ01000339.1 
             2              1              1              1              2 
VZSJ01000340.1 VZSJ01000341.1 VZSJ01000343.1 VZSJ01000344.1 VZSJ01000345.1 
             1              1              1              5              1 
VZSJ01000346.1 VZSJ01000347.1 VZSJ01000352.1 VZSJ01000354.1 VZSJ01000370.1 
             2              3              2              1              4 
VZSJ01000374.1 VZSJ01000377.1 VZSJ01000378.1 VZSJ01000380.1 VZSJ01000381.1 
             2              4              1              4              1 
VZSJ01000382.1 VZSJ01000388.1 VZSJ01000393.1 VZSJ01000395.1 VZSJ01000397.1 
             4              1              3              1              1 
VZSJ01000400.1 VZSJ01000407.1 VZSJ01000409.1 VZSJ01000415.1 VZSJ01000419.1 
             2              3              1              1              1 
VZSJ01000420.1 VZSJ01000448.1 VZSJ01000449.1 VZSJ01000453.1 VZSJ01000457.1 
             1              1              2              1             24 
VZSJ01000461.1 VZSJ01000471.1 VZSJ01000473.1 VZSJ01000481.1 VZSJ01000482.1 
             1              1              2              2              1 
VZSJ01000484.1 VZSJ01000487.1 VZSJ01000488.1 VZSJ01000490.1 VZSJ01000506.1 
             2              2              2              2              1 
VZSJ01000509.1 VZSJ01000513.1 VZSJ01000514.1 VZSJ01000515.1 VZSJ01000519.1 
             2              6              1              1              3 
VZSJ01000530.1 VZSJ01000531.1 VZSJ01000534.1 VZSJ01000544.1 VZSJ01000550.1 
             2              1              1              1              2 
VZSJ01000558.1 VZSJ01000562.1 VZSJ01000571.1 VZSJ01000590.1 VZSJ01000594.1 
             1              4              2              2              1 
VZSJ01000633.1 VZSJ01000638.1 VZSJ01000681.1 VZSJ01000694.1 VZSJ01000695.1 
             1              1              1              1              1 
VZSJ01000704.1 VZSJ01000729.1 VZSJ01000736.1 VZSJ01000752.1 VZSJ01000787.1 
             1              1              4              2              1 
VZSJ01000789.1 VZSJ01000795.1 VZSJ01000870.1 VZSJ01000876.1 VZSJ01000887.1 
             2              2              4              1              1 
VZSJ01000912.1 VZSJ01000913.1 VZSJ01000918.1 VZSJ01000938.1 VZSJ01001176.1 
             1              1              1              1              1 
VZSJ01001246.1 VZSJ01001980.1 VZSJ01001998.1 VZSJ01002064.1 VZSJ01002223.1 
             1              1              2              7              1 
VZSJ01002480.1 VZSJ01002611.1 VZSJ01002658.1 VZSJ01002689.1 VZSJ01002752.1 
             1              1              1              1              1 
VZSJ01002872.1 VZSJ01003041.1 VZSJ01003048.1 VZSJ01003085.1 VZSJ01003149.1 
             1              1              2              1              1 
VZSJ01003329.1 VZSJ01003414.1 VZSJ01003480.1 VZSJ01003492.1 VZSJ01003536.1 
             2              1              1              1              2 
VZSJ01003562.1 VZSJ01003591.1 VZSJ01003615.1 VZSJ01003636.1 VZSJ01003780.1 
             1              1              2              1              2 
VZSJ01003824.1 VZSJ01003906.1 VZSJ01003995.1 VZSJ01004006.1 VZSJ01004783.1 
             3              1              1              1              1 
VZSJ01005084.1 VZSJ01005133.1 VZSJ01005363.1 VZSJ01005398.1 VZSJ01005465.1 
             1              1              2              1              1 
VZSJ01005512.1 VZSJ01005618.1 VZSJ01012143.1 VZSJ01041603.1 VZSJ01043367.1 
             1              2              2              1              2 
Code
#overwhelmingly, the fixed differences are coming from two scaffolds: "VZSJ01000457.1" (24) and "VZSJ01000270.1" (11). Let's look at the dynamics of those scaffolds:

0.2 subset the data

Code
### plot without these two over-represented scaffolds included
auto.fixed.vcf<-fixed.vcf[fixed.vcf@fix[,1] != "VZSJ01000457.1" & fixed.vcf@fix[,1] != "VZSJ01000270.1",]
auto.fixed.vcf
***** Object of Class vcfR *****
138 samples
128 CHROMs
231 variants
Object size: 7.2 Mb
3.89 percent missing data
*****        *****         *****
Code
# Calculate hybrid index and heterozygosity for each sample. Values are returned in a data.frame
hi.het.test <- hybridIndex(vcfR = auto.fixed.vcf, pm = pm, p1 = "P0", p2 = "P1")
[1] "calculating hybrid indices and heterozygosities based on 231 sites"
Code
# View triangle plot
triangulaR::triangle.plot(hi.het.test)
Warning in geom_segment(aes(x = 0.5, xend = 1, y = 1, yend = 0), color = "black"): All aesthetics have length 1, but the data has 138 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Warning in geom_segment(aes(x = 0, xend = 0.5, y = 0, yend = 1), color = "black"): All aesthetics have length 1, but the data has 138 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Warning in geom_segment(aes(x = 0, xend = 1, y = 0, yend = 0), color = "black"): All aesthetics have length 1, but the data has 138 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.

Code
### plot just those scaffolds
z.fixed.vcf<-fixed.vcf[fixed.vcf@fix[,1] == "VZSJ01000457.1" | fixed.vcf@fix[,1] == "VZSJ01000270.1",]
z.fixed.vcf
***** Object of Class vcfR *****
138 samples
2 CHROMs
35 variants
Object size: 4.9 Mb
2.795 percent missing data
*****        *****         *****
Code
# Calculate hybrid index and heterozygosity for each sample. Values are returned in a data.frame
hi.het.test <- hybridIndex(vcfR = z.fixed.vcf, pm = pm, p1 = "P0", p2 = "P1")
[1] "calculating hybrid indices and heterozygosities based on 35 sites"
Code
# View triangle plot
triangulaR::triangle.plot(hi.het.test)
Warning in geom_segment(aes(x = 0.5, xend = 1, y = 1, yend = 0), color = "black"): All aesthetics have length 1, but the data has 138 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Warning in geom_segment(aes(x = 0, xend = 0.5, y = 0, yend = 1), color = "black"): All aesthetics have length 1, but the data has 138 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Warning in geom_segment(aes(x = 0, xend = 1, y = 0, yend = 0), color = "black"): All aesthetics have length 1, but the data has 138 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.

Code
hist(hi.het.test$heterozygosity)

Code
#make genotype plot for first scaffold
test.only<-fixed.vcf[fixed.vcf@fix[,1] == "VZSJ01000457.1",]
#isolate gt matrix
gt<-as.data.frame(t(extract.gt(test.only)))
gt$sample<-rownames(gt)
dat3 <- melt(gt, id.var = 'sample')
#plot
ggplot(dat3, aes(variable, sample)) +
  geom_tile(aes(fill = value), colour = "white") +
  scale_fill_manual(values=c("lightpink", "red", "black")) +
  theme(axis.text.y = element_text(size = 5))

Code
#repeat
test.only<-fixed.vcf[fixed.vcf@fix[,1] == "VZSJ01000270.1",]
#isolate gt matrix
gt<-as.data.frame(t(extract.gt(test.only)))
gt$sample<-rownames(gt)
dat3 <- melt(gt, id.var = 'sample')
#plot
ggplot(dat3, aes(variable, sample)) +
  geom_tile(aes(fill = value), colour = "white") +
  scale_fill_manual(values=c("lightpink", "red", "black")) +
  theme(axis.text.y = element_text(size = 5))

Code
#create sample sheet with Z chrom inversion state included
zinv<-data.frame(id=colnames(test.only@gt)[-1],
                 inv.state=c("RB","BH",rep("RB",4),"BH","het","BH",rep("het",3),
                             rep("RB",13),"het",rep("RB",28),"backcrossed RB","RB","het",
                             "RB","RB",rep("BH",17),"het",rep("BH",9),"het",
                             rep("BH",46),"RB",rep("BH",4)))

The fact that all samples seem to be either completely homozygous or completely heterozygous for fixed differences on these scaffolds suggests that this is possibly an inversion, where the whole region is inherited as a haplotype with no recombination. These two plots are highly similar, suggesting that these two scaffolds are potentially linked. To assess this, I am going BLAST them to the chromosome assembled Zebra Finch reference genome.

0.3 blast results

Code
#Blasting > 100k base-pairs of each of these two interesting scaffolds against the Zebra Finch reference genome revealed highly significant hits for both scaffolds on the Z chromosome:
knitr::include_graphics(c("/Users/devonderaad/Desktop/blast1.png"))

Code
knitr::include_graphics(c("/Users/devonderaad/Desktop/blast2.png"))

Code
#the scaffolds span at least 4 Mb, from 59 Mb to 63 Mb on the Zebra finch Z scaffold, suggesting a large inversion present in the center of the Z chromosome for these two Pheucticus taxa

0.4 make trees from the putative inverted region

Code
#Isolate Z
z.only<-v[v@fix[,1] == "VZSJ01000457.1" | v@fix[,1] == "VZSJ01000270.1",]
#get info for this dataset
z.only
***** Object of Class vcfR *****
138 samples
2 CHROMs
554 variants
Object size: 10.2 Mb
2.814 percent missing data
*****        *****         *****
Code
#z.only<-min_mac(z.only,2)
#convert to genlight
gen<-vcfR2genlight(z.only)
#fix sample names to fit in <= 10 characters
gen@ind.names
  [1] "P_hybrid_44696"         "P_hybrid_44703"         "P_hybrid_44707"        
  [4] "P_hybrid_44708"         "P_hybrid_44709"         "P_hybrid_44712"        
  [7] "P_hybrid_44762"         "P_hybrid_44771"         "P_hybrid_44781"        
 [10] "P_hybrid_45171"         "P_hybrid_45173"         "P_hybrid_45174"        
 [13] "P_ludovicianus_11998"   "P_ludovicianus_21721"   "P_ludovicianus_25286"  
 [16] "P_ludovicianus_26595"   "P_ludovicianus_33988"   "P_ludovicianus_34776"  
 [19] "P_ludovicianus_34779"   "P_ludovicianus_34782"   "P_ludovicianus_34830"  
 [22] "P_ludovicianus_44704"   "P_ludovicianus_44705"   "P_ludovicianus_44706"  
 [25] "P_ludovicianus_44710"   "P_ludovicianus_44711"   "P_ludovicianus_44713"  
 [28] "P_ludovicianus_44714"   "P_ludovicianus_44715"   "P_ludovicianus_44716"  
 [31] "P_ludovicianus_44719"   "P_ludovicianus_44720"   "P_ludovicianus_44722"  
 [34] "P_ludovicianus_44723"   "P_ludovicianus_44726"   "P_ludovicianus_44727"  
 [37] "P_ludovicianus_44729"   "P_ludovicianus_44730"   "P_ludovicianus_44731"  
 [40] "P_ludovicianus_44732"   "P_ludovicianus_44733"   "P_ludovicianus_44734"  
 [43] "P_ludovicianus_44737"   "P_ludovicianus_44738"   "P_ludovicianus_44739"  
 [46] "P_ludovicianus_44740"   "P_ludovicianus_44741"   "P_ludovicianus_44742"  
 [49] "P_ludovicianus_44743"   "P_ludovicianus_44744"   "P_ludovicianus_44745"  
 [52] "P_ludovicianus_44746"   "P_ludovicianus_44747"   "P_ludovicianus_44748"  
 [55] "P_ludovicianus_44749"   "P_ludovicianus_44753"   "P_ludovicianus_44754"  
 [58] "P_ludovicianus_44761"   "P_ludovicianus_44775"   "P_melanocephalus_34890"
 [61] "P_melanocephalus_43110" "P_melanocephalus_43276" "P_melanocephalus_44346"
 [64] "P_melanocephalus_44347" "P_melanocephalus_44471" "P_melanocephalus_44651"
 [67] "P_melanocephalus_44660" "P_melanocephalus_44666" "P_melanocephalus_44676"
 [70] "P_melanocephalus_44677" "P_melanocephalus_44678" "P_melanocephalus_44679"
 [73] "P_melanocephalus_44680" "P_melanocephalus_44681" "P_melanocephalus_44683"
 [76] "P_melanocephalus_44684" "P_melanocephalus_44685" "P_melanocephalus_44686"
 [79] "P_melanocephalus_44687" "P_melanocephalus_44688" "P_melanocephalus_44689"
 [82] "P_melanocephalus_44692" "P_melanocephalus_44693" "P_melanocephalus_44694"
 [85] "P_melanocephalus_44695" "P_melanocephalus_44697" "P_melanocephalus_44699"
 [88] "P_melanocephalus_44752" "P_melanocephalus_44760" "P_melanocephalus_44763"
 [91] "P_melanocephalus_44764" "P_melanocephalus_44765" "P_melanocephalus_44766"
 [94] "P_melanocephalus_44769" "P_melanocephalus_44770" "P_melanocephalus_44772"
 [97] "P_melanocephalus_44773" "P_melanocephalus_44774" "P_melanocephalus_44776"
[100] "P_melanocephalus_44777" "P_melanocephalus_44778" "P_melanocephalus_44779"
[103] "P_melanocephalus_44780" "P_melanocephalus_44782" "P_melanocephalus_44787"
[106] "P_melanocephalus_44788" "P_melanocephalus_44789" "P_melanocephalus_44790"
[109] "P_melanocephalus_44791" "P_melanocephalus_44792" "P_melanocephalus_44793"
[112] "P_melanocephalus_44794" "P_melanocephalus_44795" "P_melanocephalus_44798"
[115] "P_melanocephalus_44799" "P_melanocephalus_44801" "P_melanocephalus_44802"
[118] "P_melanocephalus_44803" "P_melanocephalus_44804" "P_melanocephalus_44805"
[121] "P_melanocephalus_44806" "P_melanocephalus_44808" "P_melanocephalus_44809"
[124] "P_melanocephalus_44810" "P_melanocephalus_44837" "P_melanocephalus_44840"
[127] "P_melanocephalus_44843" "P_melanocephalus_44845" "P_melanocephalus_44846"
[130] "P_melanocephalus_44847" "P_melanocephalus_44848" "P_melanocephalus_44853"
[133] "P_melanocephalus_44854" "P_melanocephalus_45175" "P_melanocephalus_45200"
[136] "P_melanocephalus_45232" "P_melanocephalus_45709" "P_melanocephalus_45926"
Code
gen@ind.names<-gsub("P_hybrid_","hyb", gen@ind.names)
gen@ind.names<-gsub("P_ludovicianus_","lud", gen@ind.names)
gen@ind.names<-gsub("P_melanocephalus_","mel", gen@ind.names)
gen@ind.names
  [1] "hyb44696" "hyb44703" "hyb44707" "hyb44708" "hyb44709" "hyb44712"
  [7] "hyb44762" "hyb44771" "hyb44781" "hyb45171" "hyb45173" "hyb45174"
 [13] "lud11998" "lud21721" "lud25286" "lud26595" "lud33988" "lud34776"
 [19] "lud34779" "lud34782" "lud34830" "lud44704" "lud44705" "lud44706"
 [25] "lud44710" "lud44711" "lud44713" "lud44714" "lud44715" "lud44716"
 [31] "lud44719" "lud44720" "lud44722" "lud44723" "lud44726" "lud44727"
 [37] "lud44729" "lud44730" "lud44731" "lud44732" "lud44733" "lud44734"
 [43] "lud44737" "lud44738" "lud44739" "lud44740" "lud44741" "lud44742"
 [49] "lud44743" "lud44744" "lud44745" "lud44746" "lud44747" "lud44748"
 [55] "lud44749" "lud44753" "lud44754" "lud44761" "lud44775" "mel34890"
 [61] "mel43110" "mel43276" "mel44346" "mel44347" "mel44471" "mel44651"
 [67] "mel44660" "mel44666" "mel44676" "mel44677" "mel44678" "mel44679"
 [73] "mel44680" "mel44681" "mel44683" "mel44684" "mel44685" "mel44686"
 [79] "mel44687" "mel44688" "mel44689" "mel44692" "mel44693" "mel44694"
 [85] "mel44695" "mel44697" "mel44699" "mel44752" "mel44760" "mel44763"
 [91] "mel44764" "mel44765" "mel44766" "mel44769" "mel44770" "mel44772"
 [97] "mel44773" "mel44774" "mel44776" "mel44777" "mel44778" "mel44779"
[103] "mel44780" "mel44782" "mel44787" "mel44788" "mel44789" "mel44790"
[109] "mel44791" "mel44792" "mel44793" "mel44794" "mel44795" "mel44798"
[115] "mel44799" "mel44801" "mel44802" "mel44803" "mel44804" "mel44805"
[121] "mel44806" "mel44808" "mel44809" "mel44810" "mel44837" "mel44840"
[127] "mel44843" "mel44845" "mel44846" "mel44847" "mel44848" "mel44853"
[133] "mel44854" "mel45175" "mel45200" "mel45232" "mel45709" "mel45926"
Code
pop(gen)<-gen@ind.names
#assign populations (a StaMPP requirement)
gen@pop<-as.factor(gen@ind.names)
#generate pairwise divergence matrix
sample.div <- stamppNeisD(gen, pop = FALSE)
#export for splitstree
#stamppPhylip(distance.mat=sample.div, file="~/Desktop/grosbeak.rad/z.only.splits.txt")
library(poppr)
This is poppr version 2.9.6. To get started, type package?poppr
OMP parallel support: unavailable
Code
ploidy(gen)<-2
genetic_distance_matrix <- poppr::bitwise.dist(gen, mat = TRUE)
#plot the two approaches to calculating genetic divergence
plot(nj(genetic_distance_matrix), type="unrooted")

Code
plot(nj(sample.div), type="unrooted")

Code
#both trees show the expected pattern (three discrete clusters)

0.5 add genomic ancestry info to dataframe

Code
#read in admixture results
#setwd to the admixture directory you brought in from the cluster
setwd("~/Desktop/grosbeak.rad/admixture.mac")

#read in input file
sampling<-read.table("binary_fileset.fam")[,1]
#get list of input samples in order they appear
sampling
  [1] "P_hybrid_44696"         "P_hybrid_44703"         "P_hybrid_44707"        
  [4] "P_hybrid_44708"         "P_hybrid_44709"         "P_hybrid_44712"        
  [7] "P_hybrid_44762"         "P_hybrid_44771"         "P_hybrid_44781"        
 [10] "P_hybrid_45171"         "P_hybrid_45173"         "P_hybrid_45174"        
 [13] "P_ludovicianus_11998"   "P_ludovicianus_21721"   "P_ludovicianus_25286"  
 [16] "P_ludovicianus_26595"   "P_ludovicianus_33988"   "P_ludovicianus_34776"  
 [19] "P_ludovicianus_34779"   "P_ludovicianus_34782"   "P_ludovicianus_34830"  
 [22] "P_ludovicianus_44704"   "P_ludovicianus_44705"   "P_ludovicianus_44706"  
 [25] "P_ludovicianus_44710"   "P_ludovicianus_44711"   "P_ludovicianus_44713"  
 [28] "P_ludovicianus_44714"   "P_ludovicianus_44715"   "P_ludovicianus_44716"  
 [31] "P_ludovicianus_44719"   "P_ludovicianus_44720"   "P_ludovicianus_44722"  
 [34] "P_ludovicianus_44723"   "P_ludovicianus_44726"   "P_ludovicianus_44727"  
 [37] "P_ludovicianus_44729"   "P_ludovicianus_44730"   "P_ludovicianus_44731"  
 [40] "P_ludovicianus_44732"   "P_ludovicianus_44733"   "P_ludovicianus_44734"  
 [43] "P_ludovicianus_44737"   "P_ludovicianus_44738"   "P_ludovicianus_44739"  
 [46] "P_ludovicianus_44740"   "P_ludovicianus_44741"   "P_ludovicianus_44742"  
 [49] "P_ludovicianus_44743"   "P_ludovicianus_44744"   "P_ludovicianus_44745"  
 [52] "P_ludovicianus_44746"   "P_ludovicianus_44747"   "P_ludovicianus_44748"  
 [55] "P_ludovicianus_44749"   "P_ludovicianus_44753"   "P_ludovicianus_44754"  
 [58] "P_ludovicianus_44761"   "P_ludovicianus_44775"   "P_melanocephalus_34890"
 [61] "P_melanocephalus_43110" "P_melanocephalus_43276" "P_melanocephalus_44346"
 [64] "P_melanocephalus_44347" "P_melanocephalus_44471" "P_melanocephalus_44651"
 [67] "P_melanocephalus_44660" "P_melanocephalus_44666" "P_melanocephalus_44676"
 [70] "P_melanocephalus_44677" "P_melanocephalus_44678" "P_melanocephalus_44679"
 [73] "P_melanocephalus_44680" "P_melanocephalus_44681" "P_melanocephalus_44683"
 [76] "P_melanocephalus_44684" "P_melanocephalus_44685" "P_melanocephalus_44686"
 [79] "P_melanocephalus_44687" "P_melanocephalus_44688" "P_melanocephalus_44689"
 [82] "P_melanocephalus_44692" "P_melanocephalus_44693" "P_melanocephalus_44694"
 [85] "P_melanocephalus_44695" "P_melanocephalus_44697" "P_melanocephalus_44699"
 [88] "P_melanocephalus_44752" "P_melanocephalus_44760" "P_melanocephalus_44763"
 [91] "P_melanocephalus_44764" "P_melanocephalus_44765" "P_melanocephalus_44766"
 [94] "P_melanocephalus_44769" "P_melanocephalus_44770" "P_melanocephalus_44772"
 [97] "P_melanocephalus_44773" "P_melanocephalus_44774" "P_melanocephalus_44776"
[100] "P_melanocephalus_44777" "P_melanocephalus_44778" "P_melanocephalus_44779"
[103] "P_melanocephalus_44780" "P_melanocephalus_44782" "P_melanocephalus_44787"
[106] "P_melanocephalus_44788" "P_melanocephalus_44789" "P_melanocephalus_44790"
[109] "P_melanocephalus_44791" "P_melanocephalus_44792" "P_melanocephalus_44793"
[112] "P_melanocephalus_44794" "P_melanocephalus_44795" "P_melanocephalus_44798"
[115] "P_melanocephalus_44799" "P_melanocephalus_44801" "P_melanocephalus_44802"
[118] "P_melanocephalus_44803" "P_melanocephalus_44804" "P_melanocephalus_44805"
[121] "P_melanocephalus_44806" "P_melanocephalus_44808" "P_melanocephalus_44809"
[124] "P_melanocephalus_44810" "P_melanocephalus_44837" "P_melanocephalus_44840"
[127] "P_melanocephalus_44843" "P_melanocephalus_44845" "P_melanocephalus_44846"
[130] "P_melanocephalus_44847" "P_melanocephalus_44848" "P_melanocephalus_44853"
[133] "P_melanocephalus_44854" "P_melanocephalus_45175" "P_melanocephalus_45200"
[136] "P_melanocephalus_45232" "P_melanocephalus_45709" "P_melanocephalus_45926"
Code
#read in all ten runs and save each dataframe in a list
runs<-list()
#read in log files
for (i in 1:10){
  runs[[i]]<-read.table(paste0("binary_fileset.", i, ".Q"))
}

#isolate run 2 (best according to CV)
run2<-runs[[2]]
#add sample info in the correct order (same as input vcf)
run2$sample<-colnames(v@gt)[-1]

#reorder
samps$sample_id == run2$sample #check if sample info table order matches the vcf
  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [37] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [49] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [85] FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
 [97] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[133] FALSE FALSE FALSE FALSE FALSE FALSE
Code
samps<-samps[match(run2$sample,samps$sample_id),] #use 'match' to match orders
samps$sample_id == run2$sample #check if this worked
  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[106] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[136] TRUE TRUE TRUE
Code
#add q-values to the sampling data.frame now that orders match
samps$mac.lud.q<-run2$V1
samps$mac.mel.q<-run2$V2

#check to see that order matches with the z-inversion dataframe I made above
samps$sample_id == zinv$id #check if this worked
  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[106] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[136] TRUE TRUE TRUE
Code
#if all above are true, add in z-inversion haplotype to the dataframe
samps$zinv.hap<-zinv$inv.state

0.6 assess inversion dynamics

Code
#how does genomic ancestry correspond to inversion haplotype in females
hist(samps$mac.mel.q[samps$zinv.hap == "BH" & samps$sex == "female"], breaks=50)

Code
hist(samps$mac.mel.q[samps$zinv.hap == "RB" & samps$sex == "female"], breaks=50)

Code
#what about in males?
hist(samps$mac.mel.q[samps$zinv.hap == "BH" & samps$sex == "male"], breaks=50)

Code
hist(samps$mac.mel.q[samps$zinv.hap == "RB" & samps$sex == "male"], breaks=50)

Code
hist(samps$mac.mel.q[samps$zinv.hap == "het"], breaks=50)

Code
#how does mitochondrial ancestry correspond to inversion haplotype in females
samps$mtDNA[samps$zinv.hap == "BH" & samps$sex == "female"]
[1] NA  1  1 NA
Code
samps$mtDNA[samps$zinv.hap == "RB" & samps$sex == "female"]
[1] NA NA NA NA  0  0  0
Code
#how does mitochondrial ancestry correspond to inversion haplotype in females
samps$mtDNA[samps$zinv.hap == "BH" & samps$sex == "male"]
 [1]  1  1  1 NA NA  1  1  1  1  1  1  1  1  1  1  1  1 NA  1  1  1  1  1  1  1
[26]  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
[51]  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 NA NA
Code
samps$mtDNA[samps$zinv.hap == "RB" & samps$sex == "male"]
 [1]  1  0  0  0  0 NA NA NA NA NA  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
[26]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
Code
samps$mtDNA[samps$zinv.hap == "het" & samps$sex == "male"]
[1] 0 0 1 1 0 0 1 0
Code
#how does genomic ancestry correspond to mitochondrial haplotype in females
hist(samps$mac.mel.q[samps$mtDNA == 0 & samps$sex == "female"], breaks=50)

Code
hist(samps$mac.mel.q[samps$mtDNA == 1 & samps$sex == "female"], breaks=50)

Code
#how does genomic ancestry correspond to mitochondrial haplotype in males
hist(samps$mac.mel.q[samps$mtDNA == 0 & samps$sex == "male"], breaks=50)

Code
hist(samps$mac.mel.q[samps$mtDNA == 1 & samps$sex == "male"], breaks=50)

Code
#check out the one weird sample with 50/50 ancestry, 50/50 phenotype, RB Z inversion and BH mtDNA
samps[samps$zinv.hap == "RB" & samps$mtDNA == 1,]
          sample_id institution catalognumber field.number  sex ossification
42   P_hybrid_44696        DMNS         44696     GMS-2277 male          100
NA             <NA>        <NA>            NA         <NA> <NA>         <NA>
NA.1           <NA>        <NA>            NA         <NA> <NA>         <NA>
NA.2           <NA>        <NA>            NA         <NA> <NA>         <NA>
NA.3           <NA>        <NA>            NA         <NA> <NA>         <NA>
NA.4           <NA>        <NA>            NA         <NA> <NA>         <NA>
NA.5           <NA>        <NA>            NA         <NA> <NA>         <NA>
NA.6           <NA>        <NA>            NA         <NA> <NA>         <NA>
NA.7           <NA>        <NA>            NA         <NA> <NA>         <NA>
NA.8           <NA>        <NA>            NA         <NA> <NA>         <NA>
     notes verbatimeventdate stateprovince                verbatimlocality
42                     39230  South Dakota SD; Lyman County; Oacoma 4mi NE
NA    <NA>                NA          <NA>                            <NA>
NA.1  <NA>                NA          <NA>                            <NA>
NA.2  <NA>                NA          <NA>                            <NA>
NA.3  <NA>                NA          <NA>                            <NA>
NA.4  <NA>                NA          <NA>                            <NA>
NA.5  <NA>                NA          <NA>                            <NA>
NA.6  <NA>                NA          <NA>                            <NA>
NA.7  <NA>                NA          <NA>                            <NA>
NA.8  <NA>                NA          <NA>                            <NA>
     decimallatitude decimallongitude
42           43.7264        -99.47388
NA                NA               NA
NA.1              NA               NA
NA.2              NA               NA
NA.3              NA               NA
NA.4              NA               NA
NA.5              NA               NA
NA.6              NA               NA
NA.7              NA               NA
NA.8              NA               NA
                                          scientificname site recordnumber
42   Pheucticus melanocephalus x Pheucticus ludovicianus    7      GMS2277
NA                                                  <NA>   NA         <NA>
NA.1                                                <NA>   NA         <NA>
NA.2                                                <NA>   NA         <NA>
NA.3                                                <NA>   NA         <NA>
NA.4                                                <NA>   NA         <NA>
NA.5                                                <NA>   NA         <NA>
NA.6                                                <NA>   NA         <NA>
NA.7                                                <NA>   NA         <NA>
NA.8                                                <NA>   NA         <NA>
     day.year weight Ltestislength Ltestiswidth Rtestislength Rtestiswidth
42        148   44.3             6            5             9            6
NA         NA     NA            NA           NA            NA           NA
NA.1       NA     NA            NA           NA            NA           NA
NA.2       NA     NA            NA           NA            NA           NA
NA.3       NA     NA            NA           NA            NA           NA
NA.4       NA     NA            NA           NA            NA           NA
NA.5       NA     NA            NA           NA            NA           NA
NA.6       NA     NA            NA           NA            NA           NA
NA.7       NA     NA            NA           NA            NA           NA
NA.8       NA     NA            NA           NA            NA           NA
     total.testis.area mtDNA bill.width bill.length tarsus wing nape back rump
42                  84     1      10.48       13.25  19.97   99    1    0    1
NA                  NA    NA         NA          NA     NA   NA   NA   NA   NA
NA.1                NA    NA         NA          NA     NA   NA   NA   NA   NA
NA.2                NA    NA         NA          NA     NA   NA   NA   NA   NA
NA.3                NA    NA         NA          NA     NA   NA   NA   NA   NA
NA.4                NA    NA         NA          NA     NA   NA   NA   NA   NA
NA.5                NA    NA         NA          NA     NA   NA   NA   NA   NA
NA.6                NA    NA         NA          NA     NA   NA   NA   NA   NA
NA.7                NA    NA         NA          NA     NA   NA   NA   NA   NA
NA.8                NA    NA         NA          NA     NA   NA   NA   NA   NA
     flanks breast.underwing.coverts male.total throat.breast.side.of.neck
42        2                        2          6                         NA
NA       NA                       NA         NA                         NA
NA.1     NA                       NA         NA                         NA
NA.2     NA                       NA         NA                         NA
NA.3     NA                       NA         NA                         NA
NA.4     NA                       NA         NA                         NA
NA.5     NA                       NA         NA                         NA
NA.6     NA                       NA         NA                         NA
NA.7     NA                       NA         NA                         NA
NA.8     NA                       NA         NA                         NA
     yellow streaking female.total distance passed.genomic.filtering mac.lud.q
42       NA        NA           NA      378                     TRUE  0.499079
NA       NA        NA           NA       NA                       NA        NA
NA.1     NA        NA           NA       NA                       NA        NA
NA.2     NA        NA           NA       NA                       NA        NA
NA.3     NA        NA           NA       NA                       NA        NA
NA.4     NA        NA           NA       NA                       NA        NA
NA.5     NA        NA           NA       NA                       NA        NA
NA.6     NA        NA           NA       NA                       NA        NA
NA.7     NA        NA           NA       NA                       NA        NA
NA.8     NA        NA           NA       NA                       NA        NA
     mac.mel.q zinv.hap
42    0.500921       RB
NA          NA     <NA>
NA.1        NA     <NA>
NA.2        NA     <NA>
NA.3        NA     <NA>
NA.4        NA     <NA>
NA.5        NA     <NA>
NA.6        NA     <NA>
NA.7        NA     <NA>
NA.8        NA     <NA>
Code
#check out the samples with het z inversions
samps[samps$zinv.hap == "het",]
                 sample_id institution catalognumber field.number  sex
107         P_hybrid_44771        DMNS         44771     GMS-2356 male
93          P_hybrid_45171        DMNS         45171    GMS-2341B male
94          P_hybrid_45173        DMNS         45173    GMS-2342B male
96          P_hybrid_45174        DMNS         45174   GMS-2344_1 male
55    P_ludovicianus_44711        DMNS         44711     GMS-2292 male
92    P_ludovicianus_44754        DMNS         44754     GMS-2335 male
33  P_melanocephalus_44685        DMNS         44685     GMS-2266 male
44  P_melanocephalus_44699        DMNS         44699     GMS-2280 male
    ossification notes verbatimeventdate stateprovince
107          100                   39238  South Dakota
93           100                   39237  South Dakota
94           100                   39237  South Dakota
96           100                   39237  South Dakota
55           100                   39234  South Dakota
92           100                   39236  South Dakota
33           100                   39229  South Dakota
44           100                   39231  South Dakota
                               verbatimlocality decimallatitude
107              SD; Lyman County; Oacoma 6mi N        43.70914
93   SD; Charles Mix County; Pickstown 20mi ESE        43.16586
94            SD; Gregory County; Platte 20mi E        43.37610
96            SD; Gregory County; Platte 20mi E        43.37347
55  SD; Charles Mix County; Greenwood 1.5mi ESE        42.96253
92   SD; Charles Mix County; Pickstown 20mi ESE        43.16900
33        SD; Mellette County; Belvidere 4mi NW        43.79035
44               SD; Lyman County; Oacoma 9mi N        43.65264
    decimallongitude                                      scientificname site
107        -99.43407 Pheucticus melanocephalus x Pheucticus ludovicianus    7
93         -98.81457 Pheucticus melanocephalus x Pheucticus ludovicianus    9
94         -99.16713 Pheucticus melanocephalus x Pheucticus ludovicianus    8
96         -99.16558 Pheucticus melanocephalus x Pheucticus ludovicianus    8
55         -98.46502                             Pheucticus ludovicianus   10
92         -98.80748                             Pheucticus ludovicianus    9
33        -101.21631                           Pheucticus melanocephalus    3
44         -99.46862                           Pheucticus melanocephalus    7
    recordnumber day.year weight Ltestislength Ltestiswidth Rtestislength
107      GMS2356      156   49.4            11            8            10
93      GMS2341B      155   44.8            12           10            10
94      GMS2342B      155   51.5            14            9            11
96     GMS2344_1      155   45.2            12            8             8
55       GMS2292      152   40.3            10            8             9
92       GMS2335      154   48.9            11            7             8
33       GMS2266      147   48.6            11           10            10
44       GMS2280      149   47.4            14            8            12
    Rtestiswidth total.testis.area mtDNA bill.width bill.length tarsus wing
107            8               168     0      11.56       13.49  22.36  104
93             9               210     0      10.33       12.90  19.96  101
94             9               225     1      10.48       13.66  20.25  102
96             8               160     1      10.83       12.32  21.08   99
55             7               143     0       9.30       12.07  19.35   96
92             7               133     0      10.13       12.88  21.49   96
33             9               200     1      10.59       13.62  19.62  108
44             9               220     0      10.31       13.55  19.33  104
    nape back rump flanks breast.underwing.coverts male.total
107    1    1    2      1                        1          6
93     1    1    1      1                        1          5
94     2    2    2      2                        3         11
96     0    0    1      1                        4          6
55     0    0    0      0                        0          0
92     1    1    1      0                        0          3
33     2    2    2      2                        4         12
44     2    2    2      1                        3         10
    throat.breast.side.of.neck yellow streaking female.total distance
107                         NA     NA        NA           NA      378
93                          NA     NA        NA           NA      460
94                          NA     NA        NA           NA      425
96                          NA     NA        NA           NA      425
55                          NA     NA        NA           NA      500
92                          NA     NA        NA           NA      460
33                          NA     NA        NA           NA      230
44                          NA     NA        NA           NA      378
    passed.genomic.filtering mac.lud.q mac.mel.q zinv.hap
107                     TRUE  0.491159  0.508841      het
93                      TRUE  0.718746  0.281254      het
94                      TRUE  0.203175  0.796825      het
96                      TRUE  0.127704  0.872296      het
55                      TRUE  0.892315  0.107685      het
92                      TRUE  0.928835  0.071165      het
33                      TRUE  0.182445  0.817555      het
44                      TRUE  0.372828  0.627172      het
Code
#check out the female samples
samps[samps$sex == "female",]
                 sample_id institution catalognumber field.number    sex
2     P_ludovicianus_21721          KU        114317              female
3     P_ludovicianus_25286          KU        117253              female
7     P_ludovicianus_34779          KU        134038              female
8     P_ludovicianus_34782          KU        134037              female
58    P_ludovicianus_44714        DMNS         44714     GMS-2295 female
59    P_ludovicianus_44715        DMNS         44715     GMS-2296 female
71    P_ludovicianus_44730        DMNS         44730     GMS-2311 female
89    P_ludovicianus_44749        DMNS         44749     GMS-2330 female
10  P_melanocephalus_34890        DMNS         34890              female
28  P_melanocephalus_44679        DMNS         44679     GMS-2260 female
112 P_melanocephalus_44776        DMNS         44776     GMS-2361 female
156 P_melanocephalus_45232        DMNS         45232     GMS-2810 female
    ossification             notes verbatimeventdate stateprovince
2            100 no phenotype info             39952        Kansas
3            100 no phenotype info             39944        Kansas
7            100 no phenotype info             43233        Kansas
8            100 no phenotype info             43222        Kansas
58           100                               39234  South Dakota
59           100                               39234  South Dakota
71           100                               39235  South Dakota
89           100                               39236  South Dakota
10             0 no phenotype info             41488      Colorado
28           100                               39221  South Dakota
112          100                               39239  South Dakota
156          100 no phenotype info             41051  South Dakota
                                          verbatimlocality decimallatitude
2                          Kansas City; 5009 Rockhill Road        39.03528
3                          Kansas City; 5009 Rockhill Road        39.03528
7                         Johnson County Community College        38.92400
8                         Johnson County Community College        38.92400
58             SD; Charles Mix County; Greenwood 1.5mi ESE        42.95152
59             SD; Charles Mix County; Greenwood 1.5mi ESE        42.93722
71                       SD; Clay County; Vermillion 5mi N        42.74242
89          SD; Lincoln County; Newton Hills; Canton 3mi N        43.22350
10                         9327 Ute Drive; Golden Co 80403        39.86550
28  SD; Lawrence County; Higgins Gulch; Spearfish 1.5mi NE        44.48330
112                      SD; Lyman County; Kennebec 12mi N        43.74631
156 SD; Lawrence County; Higgins Gulch; Spearfish 1.5mi NE        44.48330
    decimallongitude            scientificname site recordnumber day.year
2          -94.57444   Pheucticus ludovicianus   12                    NA
3          -94.57444   Pheucticus ludovicianus   12                    NA
7          -94.73500   Pheucticus ludovicianus   12                    NA
8          -94.73500   Pheucticus ludovicianus   12                    NA
58         -98.45596   Pheucticus ludovicianus   10      GMS2295      152
59         -98.43894   Pheucticus ludovicianus   10      GMS2296      152
71         -96.95320   Pheucticus ludovicianus   11      GMS2311      153
89         -96.56357   Pheucticus ludovicianus   11      GMS2330      154
10        -105.28131 Pheucticus melanocephalus    0                    NA
28        -103.93330 Pheucticus melanocephalus    1      GMS2260      139
112        -99.75561 Pheucticus melanocephalus    6      GMS2361      157
156       -103.93330                              1                    NA
    weight Ltestislength Ltestiswidth Rtestislength Rtestiswidth
2       NA            NA           NA            NA           NA
3       NA            NA           NA            NA           NA
7       NA            NA           NA            NA           NA
8       NA            NA           NA            NA           NA
58    42.5            NA           NA            NA           NA
59    47.8            NA           NA            NA           NA
71    50.0            NA           NA            NA           NA
89    49.8            NA           NA            NA           NA
10      NA            NA           NA            NA           NA
28    48.2            NA           NA            NA           NA
112   51.5            NA           NA            NA           NA
156     NA            NA           NA            NA           NA
    total.testis.area mtDNA bill.width bill.length tarsus wing nape back rump
2                  NA    NA         NA          NA     NA   NA   NA   NA   NA
3                  NA    NA         NA          NA     NA   NA   NA   NA   NA
7                  NA    NA         NA          NA     NA   NA   NA   NA   NA
8                  NA    NA         NA          NA     NA   NA   NA   NA   NA
58                  0     0      10.22       11.19  21.06  101   NA   NA   NA
59                  0     0       9.84       12.55  20.99   97   NA   NA   NA
71                  0     0      10.51       13.06  19.92   96   NA   NA   NA
89                  0     0      10.42       13.90  21.64   96   NA   NA   NA
10                 NA    NA         NA          NA     NA   NA   NA   NA   NA
28                  0     1       9.69       13.05  22.20  106   NA   NA   NA
112                 0     1      10.54       14.38  21.62  103   NA   NA   NA
156                NA    NA         NA          NA     NA   NA   NA   NA   NA
    flanks breast.underwing.coverts male.total throat.breast.side.of.neck
2       NA                       NA         NA                         NA
3       NA                       NA         NA                         NA
7       NA                       NA         NA                         NA
8       NA                       NA         NA                         NA
58      NA                       NA         NA                          0
59      NA                       NA         NA                          1
71      NA                       NA         NA                          1
89      NA                       NA         NA                          0
10      NA                       NA         NA                         NA
28      NA                       NA         NA                          3
112     NA                       NA         NA                          3
156     NA                       NA         NA                         NA
    yellow streaking female.total distance passed.genomic.filtering mac.lud.q
2       NA        NA           NA       NA                     TRUE  0.999990
3       NA        NA           NA       NA                     TRUE  0.999990
7       NA        NA           NA       NA                     TRUE  0.999990
8       NA        NA           NA       NA                     TRUE  0.999990
58       0         0            0      500                     TRUE  0.999990
59       0         1            2      500                     TRUE  0.978407
71       0         0            1      623                     TRUE  0.999990
89       1         1            2      636                     TRUE  0.983500
10      NA        NA           NA       NA                     TRUE  0.000010
28       2         3            8        0                     TRUE  0.000010
112      1         3            7      354                     TRUE  0.003458
156     NA        NA           NA       NA                     TRUE  0.000010
    mac.mel.q       zinv.hap
2    0.000010             RB
3    0.000010             RB
7    0.000010             RB
8    0.000010             RB
58   0.000010             RB
59   0.021593             RB
71   0.000010             RB
89   0.016500 backcrossed RB
10   0.999990             BH
28   0.999990             BH
112  0.996542             BH
156  0.999990             BH
Code
#is male phenotype associated with z inversion haplotype?
hist(samps$male.total[samps$zinv.hap == "RB"])

Code
hist(samps$male.total[samps$zinv.hap == "BH"])

Code
#is male phenotype associated with mtDNA haplotype
hist(samps$male.total[samps$mtDNA == 0])

Code
hist(samps$male.total[samps$mtDNA == 1])

Code
#is phenotype associated with genomic ancestry
plot(samps$mac.lud.q[samps$sex == "male"], samps$male.total[samps$sex == "male"])

Code
plot(samps$mac.lud.q[samps$sex == "female"], samps$female.total[samps$sex == "female"])

Code
#how frequent is genomic versus phenotypic intermediacy?
table(samps$male.total)

 0  1  2  3  4  5  6 10 11 12 
27  3  2  2  2  3  3  1  6 64 
Code
table(round(samps$mac.lud.q, 2))

   0 0.01 0.02 0.03 0.05 0.06 0.09 0.13 0.16 0.17 0.18  0.2 0.23 0.37 0.49  0.5 
  56    7    5    3    2    1    1    2    1    1    1    1    1    1    1    1 
0.72 0.89 0.91 0.93 0.94 0.97 0.98 0.99    1 
   1    1    1    1    1    6    3    2   37 
Code
#what proportion of samples are detectably admixed?
table(samps$mac.lud.q > 0.02 & samps$mac.lud.q < 0.98)

FALSE  TRUE 
  106    32 
Code
table(samps$mac.lud.q > 0.05 & samps$mac.lud.q < 0.95)

FALSE  TRUE 
  120    18 
Code
#what sites are those detectably admixed samples from
table(samps$site[samps$mac.lud.q > 0.02 & samps$mac.lud.q < 0.98])

 3  4  6  7  8  9 10 11 
 2  1  3  7  8  6  4  1 
Code
table(samps$site[samps$mac.lud.q > 0.05 & samps$mac.lud.q < 0.95])

 3  4  6  7  8  9 10 
 1  1  1  5  5  4  1 
Code
table(samps$site)

 0  1  2  3  4  5  6  7  8  9 10 11 12 
 8 19  3 16 10  4  8  9 10  9 10 23  9 

Both mitochondrial ancestry and Z-chromosome inversion haplotype are tightly associated with genomic ancestry. Interestingly, being heterozygous for the Z-chromosome inversion appears to be OK, and the samples that are het are spread across the transect and across the ancestry distribution. But you never see a sample with mismatched major parent ancestry and homozygous Z inversion (excluding the one 49.9/50.1 sample). Similarly, you only see one really mismatched major parent mitochondrial combination.

0.7 plot the divergence landscape

Code
#check how genotypes are encoded
table(extract.gt(v))

    0/0     0/1     1/1 
6444858  359942  102641 
Code
#extract genotype matrix
gtmat<-as.data.frame(extract.gt(v))
#recode matrix
gtmat[gtmat == "0/0"]<-0
gtmat[gtmat == "0/1"]<-1
gtmat[gtmat == "1/1"]<-2

#make sure order matches between sample sheet and genotype matrix
colnames(gtmat) == samps$sample_id
  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
 [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[106] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[136] TRUE TRUE TRUE
Code
#open up an empty vector to hold FST results
FST<-c()

#Use a for loop to calculate FST for each SNP in the genotype matrix
for(i in 1:nrow(gtmat)){
  #calc HT
    #q-bar = global derived allele frequency
    qbar<-sum(as.numeric(gtmat[i,samps$site == 0 | samps$site == 12]), na.rm = T)/(2*sum(!is.na(gtmat[i,samps$site == 0 | samps$site == 12])))
    
    #insert if statement to catch if the SNP has no variation present among the pops of interest
    if(qbar == 0 | qbar == 1){FST[i]<-"no variation"}
    else{
      
    #p-bar = global reference allele frequency
    pbar<-((2*sum(!is.na(gtmat[i,samps$site == 0 | samps$site == 12]))) - sum(as.numeric(gtmat[i,samps$site == 0 | samps$site == 12]), na.rm = T)) / (2*sum(!is.na(gtmat[i,samps$site == 0 | samps$site == 12])))
    
    #HT = 2 * p-bar * q-bar
    HT<-2*pbar*qbar

  #calc HS
    #calculate subpopulation allele frequencies 
    #p0 = subpopulation 0 alternate allele frequency
    p0<-sum(as.numeric(gtmat[i,samps$site == 0]), na.rm = T)/(2*sum(!is.na(gtmat[i,samps$site == 0])))
    #p12 = subpopulation 12 alternate allele frequency
    p12<-sum(as.numeric(gtmat[i,samps$site == 12]), na.rm = T)/(2*sum(!is.na(gtmat[i,samps$site == 12])))
    #q0 = subpopulation 0 reference allele frequency
    q0<-1-p0
    #q12 = subpopulation 12 reference allele frequency
    q12<-1-p12

    #calculate expected heterozygosities
    #Hexp1 = 1 - sum(p1^2 + q1^2)
    Hexp0<-1 - sum(p0^2 + q0^2)
    #Hexp2 = 1 - sum(p2^2 + q2^2)
    Hexp12<-1 - sum(p12^2 + q12^2)

    #HS = (Hexp1*N1 + Hexp2*N2)/Ntotal
    HS<-((Hexp0*sum(!is.na(gtmat[i,samps$site == 0]))) + (Hexp12*sum(!is.na(gtmat[i,samps$site == 12])))) / (sum(!is.na(gtmat[i,samps$site == 0 | samps$site == 12])))

  #calc and store overall FST for the given SNP
    #FST = (HT-HS)/HT
    FST[i]<-(HT-HS)/HT
    }
}

#quick double check that results make sense
table(FST == "no variation") #less than half of SNPs are actually variable (i.e., MAC > 0)

FALSE  TRUE 
23479 28116 
Code
table(FST)
FST
-1.13335267097151e-16 -2.09234339256279e-16 -3.29702595191713e-16 
                    1                     9                     1 
-3.40005801291454e-16 -9.08364292875128e-16                     0 
                    3                     1                   111 
 0.000203500203500272  0.000217013888888837  0.000238095238095498 
                    9                   199                     1 
 0.000296296296296461  0.000318877551020251  0.000341763499657948 
                    1                    30                     5 
 0.000462962962962985  0.000494071146245162  0.000496031746031218 
                   44                     2                     3 
 0.000618429189857621  0.000618429189857867  0.000686813186813569 
                    1                     1                    16 
 0.000724637681159634  0.000744047619047774  0.000754147812971423 
                    2                    13                     1 
 0.000850340136054366  0.000946969696969656  0.000987654320987844 
                    7                     1                     8 
  0.00105820105820082   0.00105820105820105   0.00105820105820166 
                   28                     2                     3 
  0.00106837606837596   0.00111607142857153    0.0012698412698411 
                    6                     7                     1 
  0.00144675925925941   0.00144927536231907   0.00160256410256374 
                    5                     3                    12 
  0.00160658055394907   0.00162337662337639    0.0017391304347832 
                    2                     5                     2 
  0.00177462289263537   0.00178253119429573   0.00178571428571434 
                    8                     2                     2 
  0.00183715461493242   0.00185185185185203   0.00189393939393954 
                   19                     1                     5 
  0.00205761316872419    0.0022321428571426    0.0022675736961455 
                    4                     6                     5 
  0.00243055555555563   0.00259740259740252    0.0026455026455027 
                    6                     2                     1 
   0.0029761904761899    0.0029761904761904   0.00305010893246185 
                    4                     1                     1 
  0.00308641975308634   0.00344827586206924   0.00347222222222223 
                    2                    23                     4 
  0.00352733686067017   0.00357142857142876   0.00366300366300374 
                    2                    10                     4 
  0.00390624999999992   0.00390624999999994   0.00392156862745098 
                    8                     5                     3 
  0.00404858299595142   0.00432900432900433   0.00444444444444416 
                    6                     3                     9 
  0.00456121145776313   0.00462962962962995   0.00468750000000002 
                   13                     1                     5 
  0.00468750000000024   0.00483091787439614   0.00496031746031741 
                    3                     6                     6 
  0.00504201680672266   0.00510204081632646   0.00510204081632647 
                    1                     2                     1 
  0.00512820512820519   0.00529100529100527   0.00534759358288742 
                    2                     7                     1 
  0.00534759358288765   0.00568181818181821   0.00571428571428571 
                    1                     1                    16 
  0.00584795321637434     0.006365740740741   0.00641025641025673 
                    2                     7                     1 
  0.00669642857142836   0.00680272108843516   0.00683760683760684 
                    4                     3                     1 
  0.00694444444444435   0.00694444444444437   0.00717703349282289 
                    1                     1                     1 
  0.00721500721500712   0.00731780167264043   0.00740740740740741 
                    4                    82                    27 
  0.00790123456790164   0.00793650793650803   0.00833333333333346 
                    1                    23                     8 
  0.00833333333333348   0.00862663906142187   0.00869565217391315 
                    1                     3                     3 
  0.00892857142857127    0.0101010101010101    0.0101750101750104 
                    2                     3                     7 
    0.010227272727273    0.0104166666666671    0.0104320987654323 
                    2                     1                     4 
   0.0105263157894738    0.0105820105820106     0.010759771629337 
                    1                    13                     3 
   0.0111408199643494    0.0112847222222222    0.0113636363636363 
                    5                    11                     2 
   0.0114468864468864    0.0114942528735632    0.0115740740740742 
                    5                    59                     1 
   0.0115900383141763    0.0116127258984403    0.0122767857142856 
                   27                     4                     4 
   0.0123456790123457     0.012475633528265    0.0124756335282652 
                    2                     1                     2 
   0.0133333333333336    0.0134205856255546    0.0136054421768705 
                   14                     3                     7 
   0.0138888888888888    0.0139204545454546    0.0142857142857142 
                    7                     1                    11 
   0.0144395078605603    0.0149342891278374    0.0151515151515152 
                    2                    83                     1 
   0.0152255639097744    0.0156249999999998              0.015625 
                    1                     2                     4 
   0.0157894736842105    0.0158371040723982    0.0158730158730158 
                    8                     2                     3 
   0.0158730158730159    0.0162037037037038    0.0166666666666667 
                    6                    14                    11 
   0.0168750000000002    0.0172532781228433     0.017543859649123 
                    1                     5                     1 
   0.0178571428571428    0.0178571428571432    0.0178571428571434 
                    1                     6                    26 
   0.0181818181818182    0.0183658933658933    0.0188297541238718 
                    4                     7                     1 
     0.01984126984127    0.0198706643151087     0.020408163265306 
                    1                     7                     1 
   0.0204081632653063    0.0204081632653064    0.0205761316872426 
                   11                     4                    17 
   0.0206851971557852     0.020833333333333    0.0208333333333332 
                    4                     3                     1 
   0.0208333333333333    0.0208333333333334    0.0217013888888885 
                   14                     2                     4 
   0.0217391304347825    0.0219587176108918    0.0220762634555742 
                    4                     7                    14 
   0.0220841959972396    0.0221608772333409    0.0222222222222222 
                    2                     1                     1 
   0.0229885057471263    0.0246924128503073     0.024805881948739 
                   96                     2                     1 
   0.0249999999999998    0.0250000000000001    0.0250896057347668 
                    2                    31                   574 
   0.0256410256410256    0.0260416666666669    0.0260869565217393 
                   15                    36                     4 
   0.0262131519274379    0.0264550264550264    0.0269360269360268 
                    2                     1                  4234 
   0.0272108843537412    0.0272222222222223    0.0277777777777777 
                    1                    14                   134 
    0.027777777777778    0.0278303872053874    0.0280092592592594 
                    6                     7                     7 
   0.0285714285714286    0.0290178571428572    0.0290516206482595 
                    5                    14                     2 
   0.0294471153846155    0.0301408179012344    0.0301408179012346 
                   10                     4                     1 
   0.0301724137931034    0.0303030303030305     0.030612244897959 
                  641                     2                     3 
   0.0310559006211177     0.032258064516129    0.0326003086419752 
                   17                  1813                     2 
   0.0327380952380948    0.0329218106995887     0.032928727046374 
                    1                     5                     1 
   0.0330687830687831    0.0340277777777778    0.0340909090909091 
                   23                     2                  2347 
    0.034285714285714    0.0342857142857146    0.0350140056022408 
                   56                     1                     1 
   0.0351562500000001    0.0352941176470588     0.035567313345091 
                    1                     7                    12 
   0.0357142857142856    0.0357142857142857    0.0361607142857142 
                    5                    43                     7 
   0.0364372469635627    0.0364583333333333    0.0366274350649349 
                    5                     1                     3 
   0.0370370370370367    0.0372670807453419    0.0375568551007151 
                  264                     3                     1 
   0.0385841836734693     0.038961038961039    0.0394088669950736 
                    3                    10                   656 
    0.039682539682539    0.0400399290150845    0.0401785714285716 
                    3                     5                    10 
   0.0412457912457914    0.0414746543778799    0.0416666666666663 
                    4                   305                     8 
   0.0417824074074073     0.042171556122449     0.042445620223398 
                   51                     2                     3 
   0.0427350427350427    0.0427350427350428     0.043097643097643 
                    1                     5                     1 
   0.0432098765432098    0.0434782608695644    0.0434782608695652 
                    1                    16                    17 
   0.0438405797101453    0.0454545454545453    0.0461538461538461 
                    1                     1                     3 
   0.0463768115942028     0.046666666666666    0.0476190476190475 
                    2                    60                     2 
   0.0476190476190477    0.0477843915343914    0.0481283422459891 
                   36                     1                     1 
   0.0482142857142856    0.0486689814814817    0.0487734487734488 
                    1                     3                     2 
   0.0492470492470493    0.0493827160493818    0.0496222527472527 
                    6                   141                    25 
   0.0506704980842913    0.0508270676691729    0.0514285714285714 
                   28                     2                    20 
   0.0517094017094019    0.0517241379310337    0.0518518518518519 
                    2                    59                   190 
   0.0520833333333334    0.0520960520960521    0.0526315789473684 
                    7                     6                     4 
   0.0527558190601671    0.0530701754385965    0.0555555555555555 
                    8                     4                     1 
   0.0555555555555557    0.0558035714285714    0.0566893424036279 
                 1408                     2                    14 
   0.0566893424036281    0.0571428571428571    0.0571428571428575 
                    1                    19                     1 
   0.0571428571428579    0.0576923076923078    0.0577887720744863 
                    2                    47                     1 
   0.0596978557504873    0.0600961538461537    0.0602083333333335 
                    3                     6                     1 
   0.0608695652173916    0.0608766233766233    0.0612244897959185 
                   13                     1                     6 
   0.0618018812463255                0.0625    0.0631313131313133 
                    5                   231                     4 
   0.0634920634920635    0.0634920634920637    0.0640000000000004 
                   10                     1                    16 
   0.0649350649350654     0.065532879818594    0.0659340659340659 
                    6                     1                     9 
   0.0666666666666667    0.0666666666666669    0.0669856459330142 
                  777                    14                     2 
   0.0679541047188105    0.0685876623376624    0.0687830687830688 
                    4                     4                     3 
   0.0688775510204081    0.0703124999999999    0.0703125000000001 
                    6                     4                   692 
   0.0714285714285713    0.0714285714285714     0.071428571428572 
                    1                     1                    26 
   0.0727272727272727    0.0736906678935667    0.0740740740740745 
                    9                     4                    15 
    0.074977817213842     0.075095785440613    0.0751537110232763 
                    8                    21                     9 
   0.0763888888888887    0.0769230769230776    0.0776502908855849 
                    2                    85                     4 
   0.0782608695652174    0.0785256410256411    0.0794477513227514 
                    8                     1                    16 
   0.0804597701149429    0.0812215724496429    0.0815217391304349 
                  100                     1                     4 
   0.0816326530612245    0.0816326530612252    0.0833333333333331 
                    1                   236                     1 
   0.0833333333333333    0.0833333333333334    0.0833829365079365 
                    9                     3                     3 
   0.0835058661145617    0.0838509316770185    0.0841750841750842 
                    4                     1                     1 
   0.0841750841750844    0.0850694444444444    0.0857142857142856 
                    2                     4                     2 
   0.0857142857142864    0.0860215053763444    0.0869565217391303 
                  113                   619                     2 
   0.0871655328798186    0.0874914559125083    0.0880208333333336 
                    4                     1                     2 
    0.089135802469136    0.0891685956790123    0.0891685956790124 
                    5                     2                     1 
   0.0892857142857144                  0.09    0.0904017857142858 
                    8                    28                     4 
   0.0909090909090911    0.0925925925925926     0.094736842105263 
                    4                     2                     3 
    0.096031746031746      0.09642094017094    0.0972222222222223 
                    6                    11                    83 
   0.0972222222222225    0.0980392156862745    0.0987755102040816 
                   24                    12                     1 
    0.098883009994121                   0.1     0.100833333333333 
                    7                     2                     1 
    0.101214574898785      0.10204081632653     0.102564102564103 
                    7                     5                    78 
    0.103448275862069     0.104184704184704     0.104637410519764 
                  328                     2                     2 
    0.105194805194805     0.106534090909091     0.107142857142857 
                    1                     8                    30 
    0.108225108225108     0.108870967741936     0.109135841836734 
                   10                   241                     2 
    0.109135841836735              0.109375     0.109750566893424 
                    2                     3                     1 
    0.111111111111111     0.111772486772487     0.111801242236025 
                   47                    21                     8 
    0.112037037037037      0.11317791005291     0.113636363636364 
                    1                    25                     5 
    0.115114795918367     0.117288961038961     0.118518518518519 
                    3                     2                   303 
    0.118822697770066                  0.12     0.120772946859903 
                    1                    33                    13 
    0.121212121212121     0.121832358674464     0.122150997150997 
                    7                     2                     4 
    0.122453703703704     0.122463768115942     0.123504273504274 
                    6                     6                     1 
                0.125     0.125744047619048     0.126482213438735 
                   17                     1                     7 
    0.126984126984127     0.127272727272727     0.127272727272728 
                  114                     3                     8 
    0.128042328042328     0.128205128205128     0.128216503992902 
                    1                     1                     4 
    0.131578947368421       0.1317738791423     0.132275132275132 
                    1                     1                     4 
    0.133004926108374     0.133282942806752     0.133333333333333 
                   36                     1                     9 
    0.133333333333334     0.134615384615385     0.135044642857143 
                    8                    57                     6 
    0.136358024691358     0.138461538461539     0.139146567717996 
                    3                     7                     1 
             0.140625     0.142857142857143     0.142908017908018 
                    3                   221                     4 
      0.1440329218107     0.146198830409357     0.146927146927147 
                   28                     1                     4 
    0.147502670940171     0.147812971342383     0.148809523809524 
                    6                     1                     3 
                 0.15     0.150554958825635     0.152173913043478 
                  133                     1                    18 
    0.153256704980843     0.155496766607878     0.155844155844156 
                  147                    16                     5 
    0.156431803490627      0.15798017771702                  0.16 
                    3                     3                    28 
    0.160079051383399     0.160714285714286      0.16304347826087 
                   11                     4                     5 
    0.163636363636364     0.166666666666667     0.166865079365079 
                   15                    39                     3 
     0.16712962962963     0.167872299382716               0.16875 
                    3                     3                     1 
    0.169388850548271     0.170068027210884     0.170556006493506 
                    2                     7                     7 
    0.170928030303031     0.173571428571429     0.173611111111111 
                    4                     1                     2 
                0.175     0.175824175824176     0.177941176470588 
                   32                    95                     1 
    0.178726035868893     0.178835978835979      0.17948717948718 
                    2                     2                    11 
    0.180515391041707     0.181481481481482     0.182449494949495 
                    1                     5                     2 
    0.183673469387755     0.184111647879764     0.184573412698413 
                   27                     2                     1 
    0.185185185185185     0.185700483091787     0.186446317657498 
                  110                     1                     6 
    0.186728395061728                0.1875      0.18796992481203 
                   10                    20                     2 
    0.188297541238718     0.188738892686261     0.190476190476191 
                    1                     4                    83 
    0.192156862745098     0.193061840120664     0.193965517241379 
                    3                     1                    62 
    0.194931773879142     0.195601851851852     0.196363636363636 
                    2                     2                     1 
    0.197530864197531     0.198380566801619                   0.2 
                    2                     8                     3 
    0.200362811791383      0.20198985042735      0.20223063973064 
                    4                    10                     4 
    0.202240896358543     0.202898550724638     0.204081632653061 
                    1                     3                     4 
    0.204545454545455     0.208695652173913     0.210069444444444 
                    2                     5                     2 
    0.212121212121212     0.215029761904762     0.215561224489796 
                   13                     2                     1 
                0.216     0.217391304347826     0.217687074829932 
                    1                    12                     4 
    0.217777777777778               0.21875     0.218762183235868 
                   10                    19                     4 
    0.221611721611722     0.222222222222222     0.223214285714286 
                    5                    29                     3 
    0.223700276331856     0.223931760204082     0.228571428571428 
                    1                     4                    41 
            0.2296875     0.230452674897119     0.230769230769231 
                   10                    80                   102 
    0.231938954765042     0.232017543859649     0.235739750445633 
                    9                     4                     2 
    0.236714975845411     0.238095238095238     0.241071428571429 
                   10                    15                    37 
    0.242126480221718     0.242424242424242     0.242712842712843 
                    1                     3                     1 
    0.243607954545454     0.246153846153846     0.249287749287749 
                    3                     3                     5 
                 0.25     0.250694444444444      0.25283950617284 
                    6                     1                     8 
    0.253968253968254      0.25568611282897     0.257142857142857 
                    5                     2                     1 
    0.259259259259259     0.262032085561497     0.264705882352941 
                   12                     1                     2 
    0.266304347826087     0.266666666666667     0.267272727272727 
                   12                     6                     1 
    0.270089285714286     0.272727272727273     0.273504273504273 
                    1                     8                    60 
    0.276315789473684     0.276734738691261     0.277777777777778 
                    1                     3                     4 
    0.278549382716049     0.279017857142857                  0.28 
                    3                     3                    63 
    0.281270200387847     0.285714285714286     0.285927854938272 
                    2                    42                     4 
    0.286549707602339     0.287423103212577     0.288089225589226 
                    1                     4                     3 
    0.288095238095238     0.289855072463768     0.290909090909091 
                    2                    20                    14 
    0.291666666666667     0.296296296296296     0.296703296703297 
                   22                     1                     6 
                  0.3     0.300662572721396     0.301247771836007 
                   11                     2                     1 
    0.301785714285714     0.304347826086957     0.305820105820106 
                    5                     8                     3 
    0.305957200694043     0.306972789115646     0.308391203703704 
                    1                     5                     6 
    0.308458737030166     0.312962962962963     0.313787737317149 
                    2                     1                     1 
     0.31578947368421            0.31640625     0.317647058823529 
                    2                     1                     2 
    0.318181818181818                  0.32     0.323978222528947 
                    8                    34                     3 
    0.325520833333333     0.326530612244898     0.327935222672065 
                    4                     1                     3 
    0.328335437710438     0.328434723171565     0.330882352941177 
                    5                     2                     1 
    0.333333333333333     0.343859649122807     0.346153846153846 
                   55                     4                    25 
    0.347222222222222     0.347556390977444     0.347826086956522 
                    1                     1                    27 
                 0.35     0.350478225478225     0.350649350649351 
                    4                     2                    10 
    0.353535353535353     0.355263157894737     0.357860331632653 
                    6                     1                     2 
                 0.36     0.360119047619048     0.360428849902534 
                    4                     1                     4 
    0.361480624638519     0.363636363636364     0.365017361111111 
                    1                    13                     2 
    0.368421052631579     0.368622448979592     0.369126043039087 
                    1                     3                     6 
     0.37037037037037                 0.375     0.377232142857143 
                   38                    15                     5 
    0.380952380952381      0.38507326007326     0.389044506691565 
                    6                     6                     1 
             0.390625     0.391304347826087      0.39360119047619 
                    2                    39                     2 
    0.395093795093795     0.396825396825397                   0.4 
                    2                     6                     7 
                0.405     0.407407407407407     0.409844054580897 
                   16                     5                     2 
    0.415384615384615     0.415584415584416     0.416666666666667 
                    1                    21                    13 
    0.417100694444444      0.41984126984127      0.42512077294686 
                    3                     4                    39 
    0.426118827160494     0.428571428571429     0.429824561403509 
                    5                     4                     3 
    0.433155080213904     0.434389140271493      0.43565867003367 
                    2                     3                     6 
               0.4375     0.442469295410472     0.444444444444444 
                    3                     4                     2 
    0.444444444444445     0.446428571428572     0.447668650793651 
                    7                     2                     4 
    0.453781512605042     0.454545454545455     0.455314422419686 
                    3                    21                     3 
    0.458333333333333     0.466666666666667     0.466709760827408 
                    1                     2                     3 
              0.46875     0.473373187658902     0.473684210526316 
                    5                     1                     8 
    0.474509803921569     0.474567901234568                  0.48 
                    5                     1                     1 
    0.484848484848485     0.485294117647059     0.489795918367347 
                   29                     4                    14 
    0.489878542510121     0.493059446000623     0.497124756335283 
                    4                     1                     5 
                  0.5     0.501736111111111     0.503105590062112 
                    7                     1                     3 
    0.504166666666667     0.506578947368421     0.508750508750509 
                    1                     8                     4 
    0.510204081632653     0.517007797270955     0.518518518518519 
                    1                     2                     1 
    0.520833333333333     0.523809523809524     0.532163742690059 
                    2                    20                     3 
    0.533333333333333     0.535714285714286      0.53804347826087 
                    4                     1                     7 
    0.542410714285714     0.549186862244898      0.55026455026455 
                    3                     2                    29 
    0.555555555555556     0.556563366087176      0.55978835978836 
                    7                     2                     2 
               0.5625     0.571428571428571     0.571428571428572 
                    7                    19                     1 
    0.576923076923077     0.583333333333333     0.583526234567901 
                    1                    11                     1 
    0.584415584415584     0.589335317460317     0.589473684210526 
                    4                     2                     4 
    0.594184027777778                   0.6     0.604938271604938 
                    2                    24                     4 
    0.613636363636364      0.61764705882353      0.62051282051282 
                    6                     2                     2 
    0.622222222222222                 0.625     0.628571428571428 
                   28                     1                     2 
    0.631578947368421       0.6400290885585     0.642857142857143 
                    6                     2                     5 
    0.647058823529412                  0.65     0.661654135338346 
                    6                     1                    14 
    0.662745098039216     0.669117647058823     0.673469387755102 
                    6                     6                     1 
    0.678557504873294     0.680555555555556      0.68421052631579 
                    4                     7                    27 
    0.686274509803922     0.694444444444445     0.696428571428571 
                    3                     1                     5 
    0.701754385964912     0.711111111111111     0.729166666666667 
                   13                     1                     2 
    0.740740740740741     0.741512345679012                  0.75 
                    3                     1                     3 
    0.750079719387755     0.761904761904762     0.762156714537667 
                    2                     7                     2 
             0.765625     0.771428571428571     0.771428571428572 
                   11                     1                     5 
    0.777777777777778      0.77782600308642                0.7875 
                   22                     6                     9 
    0.790123456790124     0.823529411764706     0.847058823529412 
                   16                     1                     3 
    0.852272727272727     0.857142857142857     0.862745098039216 
                    2                     1                     4 
    0.865384615384615     0.866666666666667     0.873949579831933 
                    4                     2                    12 
                0.875     0.879699248120301     0.881481481481481 
                    6                     1                     5 
     0.88235294117647     0.888157894736842     0.888888888888889 
                   22                     5                    34 
                    1  3.99680288865056e-16  4.87329434697457e-05 
                  266                     1                    10 
 5.66676335485757e-16  7.66812361015085e-05  8.08015513897186e-05 
                    1                     3                     3 
 9.07029478457844e-05          no variation 
                   17                 28116 
Code
#make dataframe
FST.dataframe<-data.frame(chrom=v@fix[,1],pos=v@fix[,2],FST=FST)
#isolate variable SNPs
FST.dataframe<-FST.dataframe[FST.dataframe$FST != "no variation",]
head(FST.dataframe)
           chrom    pos                FST
1 VZSJ01000159.1 151719 0.0804597701149429
2 VZSJ01000159.1 151774 0.0250896057347668
4 VZSJ01000159.1 152160 0.0229885057471263
5 VZSJ01000159.1 152159 0.0476190476190477
6 VZSJ01000159.1 152142 0.0229885057471263
7 VZSJ01000159.1 152124 0.0517241379310337
Code
#make numeric
FST.dataframe$FST<-as.numeric(FST.dataframe$FST)
head(FST.dataframe)
           chrom    pos        FST
1 VZSJ01000159.1 151719 0.08045977
2 VZSJ01000159.1 151774 0.02508961
4 VZSJ01000159.1 152160 0.02298851
5 VZSJ01000159.1 152159 0.04761905
6 VZSJ01000159.1 152142 0.02298851
7 VZSJ01000159.1 152124 0.05172414
Code
#plot histogram of FST
hist(FST.dataframe$FST)

Code
#plot divergence landscape
plot(1:nrow(FST.dataframe), FST.dataframe$FST, cex=0.1, pch=19)

Code
plot(1:nrow(FST.dataframe), FST.dataframe$FST, cex=0.1, pch=19, col=as.factor(FST.dataframe$chrom))

Code
palette(c("grey","black"))
plot(1:nrow(FST.dataframe), FST.dataframe$FST, cex=0.1, pch=19, col=as.factor(FST.dataframe$chrom))

Code
#explore
r1<-with(FST.dataframe, tapply(FST, chrom, mean))
hist(r1)

Code
r1[r1 > .5]
VZSJ01000272.1 VZSJ01000339.1 VZSJ01000345.1 VZSJ01000354.1 VZSJ01000393.1 
     0.5374647      0.6634282      1.0000000      0.5170455      0.5611313 
VZSJ01000481.1 VZSJ01000487.1 VZSJ01000490.1 VZSJ01000530.1 VZSJ01000534.1 
     0.5859914      0.6780303      1.0000000      0.7008547      1.0000000 
VZSJ01000576.1 VZSJ01000681.1 VZSJ01000736.1 VZSJ01000938.1 VZSJ01002872.1 
     0.5533772      1.0000000      0.6008742      1.0000000      0.5351562 
VZSJ01003329.1 VZSJ01003995.1 VZSJ01004006.1 VZSJ01012143.1 VZSJ01043367.1 
     0.5161290      1.0000000      0.8088235      1.0000000      0.5740285 
Code
hist(table(FST.dataframe$chrom), breaks=100)

Code
table(FST.dataframe$chrom[FST.dataframe$FST == 1])

VZSJ01000159.1 VZSJ01000241.1 VZSJ01000254.1 VZSJ01000265.1 VZSJ01000270.1 
             3              2              3              5             11 
VZSJ01000272.1 VZSJ01000279.1 VZSJ01000280.1 VZSJ01000283.1 VZSJ01000285.1 
             3              9              3              4              3 
VZSJ01000299.1 VZSJ01000301.1 VZSJ01000305.1 VZSJ01000337.1 VZSJ01000339.1 
             2              1              1              1              2 
VZSJ01000340.1 VZSJ01000341.1 VZSJ01000343.1 VZSJ01000344.1 VZSJ01000345.1 
             1              1              1              5              1 
VZSJ01000346.1 VZSJ01000347.1 VZSJ01000352.1 VZSJ01000354.1 VZSJ01000370.1 
             2              3              2              1              4 
VZSJ01000374.1 VZSJ01000377.1 VZSJ01000378.1 VZSJ01000380.1 VZSJ01000381.1 
             2              4              1              4              1 
VZSJ01000382.1 VZSJ01000388.1 VZSJ01000393.1 VZSJ01000395.1 VZSJ01000397.1 
             4              1              3              1              1 
VZSJ01000400.1 VZSJ01000407.1 VZSJ01000409.1 VZSJ01000415.1 VZSJ01000419.1 
             2              3              1              1              1 
VZSJ01000420.1 VZSJ01000448.1 VZSJ01000449.1 VZSJ01000453.1 VZSJ01000457.1 
             1              1              2              1             24 
VZSJ01000461.1 VZSJ01000471.1 VZSJ01000473.1 VZSJ01000481.1 VZSJ01000482.1 
             1              1              2              2              1 
VZSJ01000484.1 VZSJ01000487.1 VZSJ01000488.1 VZSJ01000490.1 VZSJ01000506.1 
             2              2              2              2              1 
VZSJ01000509.1 VZSJ01000513.1 VZSJ01000514.1 VZSJ01000515.1 VZSJ01000519.1 
             2              6              1              1              3 
VZSJ01000530.1 VZSJ01000531.1 VZSJ01000534.1 VZSJ01000544.1 VZSJ01000550.1 
             2              1              1              1              2 
VZSJ01000558.1 VZSJ01000562.1 VZSJ01000571.1 VZSJ01000590.1 VZSJ01000594.1 
             1              4              2              2              1 
VZSJ01000633.1 VZSJ01000638.1 VZSJ01000681.1 VZSJ01000694.1 VZSJ01000695.1 
             1              1              1              1              1 
VZSJ01000704.1 VZSJ01000729.1 VZSJ01000736.1 VZSJ01000752.1 VZSJ01000787.1 
             1              1              4              2              1 
VZSJ01000789.1 VZSJ01000795.1 VZSJ01000870.1 VZSJ01000876.1 VZSJ01000887.1 
             2              2              4              1              1 
VZSJ01000912.1 VZSJ01000913.1 VZSJ01000918.1 VZSJ01000938.1 VZSJ01001176.1 
             1              1              1              1              1 
VZSJ01001246.1 VZSJ01001980.1 VZSJ01001998.1 VZSJ01002064.1 VZSJ01002223.1 
             1              1              2              7              1 
VZSJ01002480.1 VZSJ01002611.1 VZSJ01002658.1 VZSJ01002689.1 VZSJ01002752.1 
             1              1              1              1              1 
VZSJ01002872.1 VZSJ01003041.1 VZSJ01003048.1 VZSJ01003085.1 VZSJ01003149.1 
             1              1              2              1              1 
VZSJ01003329.1 VZSJ01003414.1 VZSJ01003480.1 VZSJ01003492.1 VZSJ01003536.1 
             2              1              1              1              2 
VZSJ01003562.1 VZSJ01003591.1 VZSJ01003615.1 VZSJ01003636.1 VZSJ01003780.1 
             1              1              2              1              2 
VZSJ01003824.1 VZSJ01003906.1 VZSJ01003995.1 VZSJ01004006.1 VZSJ01004783.1 
             3              1              1              1              1 
VZSJ01005084.1 VZSJ01005133.1 VZSJ01005363.1 VZSJ01005398.1 VZSJ01005465.1 
             1              1              2              1              1 
VZSJ01005512.1 VZSJ01005618.1 VZSJ01012143.1 VZSJ01041603.1 VZSJ01043367.1 
             1              2              2              1              2 
Code
#check whether FST is elevated on the putative inversion chromosomes versus the rest of the genome
mean(FST.dataframe$FST)
[1] 0.09345442
Code
mean(FST.dataframe$FST[FST.dataframe$chrom == "VZSJ01000457.1"])
[1] 0.2522709
Code
mean(FST.dataframe$FST[FST.dataframe$chrom == "VZSJ01000270.1"])
[1] 0.2564168
Code
#FST is clearly elevated on those chroms 


#it seems like we don't really have the resolution to get a detailed look at the divergence landscape, which would probably require binning FST in sliding windows and better scaffolding. Best left for a future whole genome follow-up to assess how important the inversion is in separating hybrids and whether it is segregating within each pop at all.